Oscillations and patterns in interacting populations of two species 



Matti Peltomaki 1 , Martin Rost 2 , Mikko Alava 1 
1 Department of Applied Physics, Helsinki University of Technology, P.O. Box 1100, 02015 HUT, Espoo, Finland 
2 Bereich Theoretische Biologie, IZMB, Universitdt Bonn, 53115 Bonn, Germany 

Interacting populations often create complicated spatiotemporal behavior, and understanding it 
is a basic problem in the dynamics of spatial systems. We study the two-species case by simulations 
of a host-parasitoid model. In the case of co-existence, there are spatial patterns leading to noise- 
sustained oscillations. We introduce a new measure for the patterns, and explain the oscillations as 
a consequence of a timescale separation and noise. They are linked together with the patterns by 
letting the spreading rates depend on instantaneous population densities. Applications are discussed. 
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A fundamental aim in studying population dynamics 
is to understand species interactions. A paradigmatic, 
still interesting, case is that of two species where one 
feeds or lives from the other. These may be predators 
and their prey or parasitoids living on the expense of 
hosts. If the interaction is strong and if the parasitoid 
or predator is specialized, its growth will have a delayed 
negative feedback as its resource is diminished. In such 
cases, oscillations are an inherent feature 9. Many sur- 
face reactions also have similar dynamics [3( . 

The classical ways of looking at such systems assume 
fully stirred populations. The encounters between indi- 
viduals are assumed to be proportional to the product of 
their densities, analogously to the mass action principle 
in chemistry This assumption is also at the heart 
of the mean-field -like Lotka-Volterra equations. In gen- 
eral, however, spreading and interaction are restricted in 
space. In this case, correlated structures arise and the 
assumption about complete mixing no longer holds 

If the parasite abundance is small, any feedback effect 
is weak. Population sizes then show no oscillations, and 
the predating species is locally concentrated in a cluster- 
like arrangement. This has been theoretically studied in 
Ref. 0. With strong feedback spatiotemporal patterns 
emerge in a multitude of forms. These include disordered 
flame- likepatterns [a, 0, @| , and ripple- like spatiotempo- 
ral ones [9(. There is also a large body of work on sim- 
ilar patterning in individual-based models (e.g. (loj). in 
statistical physics 11, 12j], and in calcium concentration 
oscillations in living cells A paradigmatic pattern- 
forming system is the complex Ginzburg-Landau equa- 
tion (CGLE) [3], which exhibits spiral-like geometries. 

Voles in Northern Britain [HI], mussels in the Wad- 
den sea [3], and lemmings in Northern Europe 17], are 



good examples of empirical observations of such pattern- 
ing. These involve either predation or being predated. 
Spatial structures weaken the interactions since species 
tend to be aggregated within themselves. They also pro- 
vide the prey a refuge since around the prey there are less 
predators. Therefore, spatial inhomogeneity can stabilize 
the dynamics and promote coexistence 0, 0, . 

Here, we analyze the full spatiotemporal dynamics of 



two interacting species using instantaneous configura- 
tions and time series of the population densities. We 
introduce a measure for the level of patterning in such 
systems. When the patterns form, one observes persis- 
tent oscillations in the population densities. We show 
that the underlying dynamics follows a particular logic: 
it originates in the response of the rates to changes in 
instantaneous densities, and the emerging system proves 
different from the limit cycle in Lotka-Volterra systems, 
or recent developments where three-species models have 
been mapped 12, 2(| to CGLE. The present mechanism 
works by the interplay of oscillatory transients to a stable 
fixed point and stochasticity. The response of the interac- 
tion rates is due to spatial correlations. The mechanism 
is novel, and does not in fact need any non-linearities to 
work, which will become evident below based on simula- 
tions and effective equations describing them. 

We study a host-parasitoid model in discrete time and 
space. It is inspired by [2l|, but has a wider interaction 
range as in the incidence function models of metapopu- 
lation dynamics (2^]. The model describes annual host- 
parasitoid dynamics on a two-dimensional square lattice 
A 23] . At each time step, a site x can be either empty (in 
state e) or populated by a host without (state h) or with 
parasitoids (state p) . Transitions between the states are 
cyclic, e — > h — > p — > e, neglecting possible spontaneous 
deaths of non-parasitized hosts assumed to be rare, for 
simplicity. Although this means that hosts live forever if 
there are no parasitoids, this is not a serious restriction; 
in coexistence it boils down to assuming faster extinction 
for the parasitoids than for the hosts. The model can also 
be described as an SIR model with rebirth. 

At each site transition probabilities depend on the sur- 
rounding populations through the connectivity 

J a (x,t) = MI X - X 'D Xa{x',t) (1) 

x'eA 

of site x with respect to species a {h or p). Here 
Xa(x) = 1 if the state of x is a, and else. The ker- 
nel k a has an exponential decay with the scale w a and 
is normalized by J2 x eA ^a(l x l) = 1- Dispersal lengths 
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and 



are chosen since these are biologically motivated [25 
lead to generalizability. 

In a timestcp, the transition e — > h takes place with 
probability min(l, A/jl/i) and h — > p with probability 
min(l, Xplp) (in the parameter range of interest, \ a I a < 
1 practically always). Parasitized hosts may die (p — > e) 
with probability 5 irrespective of the surroundings. Note 
that they do not reproduce. Periodic boundary condi- 
tions and parallel updates are used. There are two ab- 
sorbing states, an empty lattice (e) and one full of hosts. 
Fig. [1] (a) shows an example with coexistence. 
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FIG. 1: (a) Patterned co-existing configuration with param- 
eters Wh = 3.0, w p — 1.5, \h = 0.63, A p = 2.5, S = 0.9, and 
L = 512. Sites in e are white, h is gray, and p black, (b) 
Dominance regions of the same state, (c) The definition of 
the coarse domains. For given h and p, the first quadrant of 
the (ph, p p )-plane is divided into three regions by three lines. 
That separating h- and p-dominated regions starts from the 
average (h,p) (the black dot), goes towards increasing ph and 
p p and is such that its continuation passes through the origin. 
The other two lines form 120-degree angles with the first one 
and each other, (d) The domain wall length ratio for different 
values of Ah and A p . Other parameters are 5 = 0.9, Wh = 3.0, 
w p — 1.5, L = 512, and a = 8. 

One finds moving regions predominantly in one of 
the three states. To quantify these patterns, we de- 
fine the dominance regions (Fig. lb) as follows. By 
smoothing one obtains continuous densities p Q (x, £) — 
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Xa (x' , t) . For each site, ph (x, t) and 



/9 p (x, t) are compared to the space-time averages h and p. 
The densities are positive, lying in the first quadrant of 
M 2 , divided into three regions shown in Fig. lc. The site 
x at time t is then defined to belong to a domain accord- 
ing to the region (e, h, or p) containing (ph{x, t), p p (x, t)). 
In essence, the regions coarse-grain on a scale a > Wh. p - 
In this regime, they are insensitive to changes in a. 

The domains are separated by walls, joining at triple 
points, vortices @. A vortex has a sign +1 (—1), if one 
encounters the domains in the order ehp following a small 
cycle around the vortex counter-clockwise (clockwise). 
Pairs of vortices of opposite signs are created and annihi- 
lated together. The domains rotate around the vortices, 



which are relatively stable. In other words, the species 
invade the appropriate neighboring domains so that the 
walls rotate around the vortices. Similar structures have 
been identified earlier in related systems (e.g. frlflljo. I n 
three dimensions, the vortices generalize to strings |6(. 

First, consider static measures such as the domain wall 
length from source to sink vortex. It has an exponential 
distribution, whose mean is drawn for different parame- 
ters in Fig. Oil. Its ratio to its counterpart in uncorrelated 
random arrangements with the same densities is shown. 
Patterns and oscillations lead to walls with (£) w 100 
lattice units (l.u.). This is more than ten times larger 
than the smoothing width a and also several times that 
in the random arrangements (Random = 35 l.u.). The 
coarse-graining gives a measure of patterns distinguish- 
ing between uncorrelated and patterned states. 

Next, turn to the spatially averaged densities h t and 
Pt- Fig. [2k shows them as a function of time in three 
cases: (i) a non-patterned state with a small parasitoid 
population, (ii) a state with patterns, and (iii) a small 
subsystem (L su b = 64) out of a large system (L = 512) 
with patterns. In the patterned systems, there is a high- 
frequency oscillation matching the angular velocity of sin- 
gle vortices and a slow variation of the amplitude. Below, 
we explain both as a consequence of a timescale separa- 
tion, connect them to the patterns, and explain why the 
oscillations do not conform to the usual limit cycles. 
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FIG. 2: (a) Upper panel: the densities in the patterned state 
for L = 512 (solid lines), and a subsystem with L — 64 
(dashed lines). See Fig. [1] for parameters. Lower panel: the 
same in the homogeneous state for L = 512 and A p = 1.3. 
The parasitoid curve has been shifted for clarity, (b) The 
phase diagram of the system. (1) Parasitoids extinct, (2) non- 
oscillatory coexistence, (3) coexistence with noise-sustained 
oscillations, (4) coexistence in a limit cycle, and (23) a tran- 
sition zone between (2) and (3). Black (red) lines denote the 
boundaries for the spatial system (MF approximation). The 
boundary between (1) and (2) coincides for the two cases. 



To build a description of the dynamics of the model 
in a novel fashion using aggregated variables, consider 
Poincare maps (Fig. [3]). For large enough systems h t +i 
and pt+i are unique functions of ht and pt up to noise. 
Also by attractor reconstruction 24j we find that the full 
system with 2L 2 degrees of freedom coarse-grains into a 
two-dimensional one. The points in the maps lie close 
to a two-dimensional surface, and for large enough L 
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(with many patterns in the system) even on the tangen- 
tial plane through the average (h,p). Based on these 
numerical observations, the dynamics is linear in h t and 
Pt- 



h t +i = a/^hht + a,h, p Pt + c h 
Pt+i = a P ,hh t + a PtP p t + c p . 



(2) 



In other words, the observations imply that even though 
the dynamics is expected to be non-linear based on the 
MF approximation, in a large system with many patterns 
the possible non-linearities self-average out. 
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FIG. 3: Left: ht+i over h t and pt (black) in a perspective 
showing the planar arrangement of points and its projection 
onto the (/it,pt)-plane (gray). Right: the same for pt+i- Pa- 
rameters are as in Fig. [T] 

It is then useful to consider the expansion around the 
average - a linear iterative map. In oscillatory cases, 
its eigenvalues form a conjugated pair pe ±z ^. These are 
associated with two timescales, the period and the decay 
rate of the amplitude. Their ratio 



v = 0/|logp| 



(3) 



tells whether dynamics is oscillatory or "just noisy". 
v ^> 1 indicates patterned systems and oscillatory dy- 
namics. Fig. 2] shows the typical behaviors of the two 
kinds of dynamics observed depending on the presence or 
absence of patterns, and whether one adds noise (as ad- 
ditional Gaussian uncorrelated noise terms on the RHS 
of Eq. |2])) to the coarse-grained dynamical system to 
mimic the finite-!/ simulations of the full spatial system. 
Note that in all cases the fixed point is attractive. 

So far we have given separately a temporal and a spa- 
tial diagnosis of the pattern dynamics. Next we link these 
together. For X a I a (x,t) small, they equal the spread- 
ing probabilities, and the dynamics can be written as 

ht+i = h t +X)xeA A/ l fc/ l (x)C e/l (x, t) - X p k p (x)C h p(x,t) 

andpt+i = (l-8)p t + \ p 2xgaM x ) C hp (x,t), where the 
influence of the connectivities is expressed by the corre- 
lation functions C a p(x,t) = X^x'eA Xa(x' , t) X/?( x + 
x',t). A corresponding non-spatial approximation is 



h t +i = h t + K(ht,Pt)(l-h t -pt)ht 
Pt+i = (1 - S)p t + p(h t ,pt)h t pt . 



p(ht,pt)h t pt 
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FIG. 4: The behavior of Eqs. ([2]) in the patterned (a, b, c) and 
in the homogeneous state (d, e, f). The coefficients a a a i and 
c a are obtained from fits of Eq. @ to simulation data. In both 
cases they lead to an oscillatory convergence to the FP (a, 
d). With noise, the FP is never reached and the cases differ: 
in the homogeneous one this results in random fluctuations 
around the FP (e); in the patterned case the noise kicks the 
system out of the FP with slow and oscillatory decay and 
leads to persistent oscillations different from limit cycles (b). 
The host density as a function of time corresponding to (b, e) 
is shown in (c, f). 



This is an approximation of the usual MF form with 
the interaction parameters n{h,p) and /i(h,p) generalized 
to be arbitrary functions of the instantaneous densities. 
They can be non-linear and they do not have to conform 
to the standard MF nor to any ad-hoc approximations 
[2o| . By an expansion of Eqs. §4$ around the fixed point 
(h,p) one arrives at Eq. |(5J) with 



ah,h 
a p ,h 



l + n — 2Kh—(n+[j,)p+dhKh(l~h—p) — dh^hp 
— (K-\-(j,)h — d p K h(l—h—p) — 9p/i hp 
[ip + dhp, hp (5) 
1 — 5 + fih + dpfi hp, 



(4) 



where n, fx, and their derivatives are evaluated at the 
fixed point. The derivatives are necessary for consistency. 
The matrix elements a a . a i and the densities h and p are 
measured from the simulations. Since k and fi are pa- 
rameters, omitting the derivatives would make Eqs. ([5]) 
overdetermined and thus unsatisfiable. By keeping them, 
there are four equations and four unknowns to be solved. 

The effect of the nonzero derivatives is best illustrated 
by a phase diagram, Fig. [21 In the spatially extended 
system, there are three qualitative phases: the extinction 
of the parasitoids, non-oscillatory and oscillatory coexis- 
tence. Except for the extinction, the boundaries are not 
sharp. Instead, there is a transition zone, defined via the 
timescale ratio (Eq. ©) as the region where 1<^<4. In 
MF, there is also a fourth phase, absent here: oscillatory 
coexistence in a limit cycle. The phase structure resem- 
bles that in earlier work on a related model with only 
nearest- neighbour spreading [26], [27, 28 1. There, as well, 
oscillatory and non-oscillatory phases are recovered, the 
latter identified as a limit cycle using the pair approx- 
imation. Based on our findings, it could also be noise- 
sustained. 

Let us now compare the explained mechanism with 



4 



recent approaches. A possibility is to make the angu- 
lar velocity of the oscillation either amplitude- or phase- 
dependent 2f| 3(J. However, Eq. (T5|) does not allow for 



either dependence. Another one is to map the population 
model 1 



20J to CGLE [14j. An unstable fixed point is 
necessary for the mapping, yielding a limit cycle. 

To conclude, we have studied spatiotemporal dynam- 
ics of a model of two kinds of interacting particles, in 
biological terms hosts and parasitoids. A large para- 
sitoid population creates patterns and noisy oscillations 
of population sizes. We have introduced a new measure 
for the patterns, and explained the noisy oscillation as a 
consequence of a time-scale separation. In other words, 
even with a limit cycle at the well-mixed limit, the spa- 
tial case has stable dynamics with long-lived oscillatory 
transients. This is due to spatial correlations making the 
spreading rates functions of the instantaneous popula- 
tion densities. Since the type of oscillation determines 
its properties (e.g. the fluctuating amplitude), which in 
turn affect vulnerability to extinction, the distinction is 
important. The connection offers a shortcut to study the 
effect of, e.g. environment: it could be related directly 
to the matrix elements in ([2]), in contrast to a full form 
of the interactions, lightening the analysis. We expect 
that the observation of patterns and oscillations arising 
from local dynamics and self-averaging - in finite sys- 
tems since the noise amplitude depends on system size - 
will find other applications beyond the biology-inspired 
model. They are not restricted to only two species or 
two-dimensional systems, since the analysis can be car- 
ried out also for more complicated cases. There is no 
restriction to cyclic dynamics either, nor to discrete-time 
systems since continuous-time ones can be handled by 
considering snapshots taken at regular intervals. Fur- 
ther examples of applications include chemical reactions 
on surfaces H, and metapopulations on disordered and 
scalefree landscapes. 
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